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1. INTRODUCTION 


The accurate modeling of oceanic tides is very important for the inteipretation of oceanic measurements 
from space borne altimeters. The tides not only introduce errors into the determination of geostrophic velocity 
but they are of considerable interest in their own right. Tidal investigations of ocean basins with realistic 
topography and coastal boundaries are still in the developing stage. 

Ocean tidal studies have a long history. Several comprehensive reviews can be found in the literature, 
i.e., Hendershott and Munk (1970), Cartwright (1977), Hendershott (1973, 1977, 1981) and Schwiderski 
(1980). 

Hendershott (1977) summarized the efforts to solve the LTE for global tides for semidiurnal and diurnal 
components up to that time. The majority of the M2 calculations agree in the prediction of both a North 
and a South Atlantic amphidrome, although they vary widely in the resolution and prediction of coastal 
tides and none is able to resolve the tides over the shallow Patagonian shelf. An anti-amphidrome is predicted 
consistently in the central Indian ocean with the coastal details showing wide variation for the different models. 
In the southeast Pacific some calculations yield a confluence of cotidal lines with low amplitude while other 
models locate an amphidrome there. No data are available to resolve the discrepancy. The northwest Pacific 
is different in every model. 

More recent efforts have been directed to incorporate the results of loading and self-attraction. Accad 
and Pekeris (1978) solve the Laplace’s tidal equations for the M2 and S2 tides in the world oceans on the 
basis of a knowledge of the tidal potential alone. Tidal dissipation was taken to be limited to the coastline. 
The main purpose of their investigation was to determine the effects of tidal self-attraction and of tidal loading. 
An iterative method was developed to evaluate these secondary effects, the resulting change is of the order 
of 10% and better agreement is obtained between the theoretical and observed tides. Parke and Hendershott 
(1980) obtained global solutions to the Laplace’s tidal equations for the M2, S2 and K1 tides. They ad- 
dressed the problem of divergence of the near-resonance modes by means of test functions which are used 
to interpolate between island data in the least squares sense. These test functions are derived by solving 
Laplace’s tidal equations with ocean loading and selfgravitation in an iterative manner. The resulting represen- 
tations of the global tide are stable over at least a ±5% variation in the mean depth of the model basin 
and they conserve mass. 

Schwiderski (1980) computed the M2 global tide by solving the Laplace’s tidal equations directly, using 
a finite difference scheme in space and time; the strictly mathematical solution was modified by means of 
an interpolation technique which incorporates over 2000 tide data collected at continental and island sta- 
tions, the interpolation is accomplished by adjusting the bottom friction coefficient and by allowing in or 
out-flow across the mathematical ocean boundary. No direct comparison of observed and computed data 
is feasible since the model incorporates essentially all known data although it is possible to evaluate the 
smoothness with which the computed tide accepts or rejects data. It was found that the interpolation tech- 
nique permits a check of the reality of both the tide model and the tidal data input. 

Platzman (1978, 1981,1984) has computed a range of normal modes for the world oceans. He used them 
to synthesize some of the diurnal and semidiurnal tides. He decomposed the transport vector by means of 
the Stokes/Helmholtz potentials but did not determine the velocity potential and stream function eigenfunc- 
tions explicitly (as in this paper), but proceeded directly to the normal mode solution of the L.T.E. The 
linearized primitive equations were discretized by means of first-order, piecewise-linear finite elements with 
an average grid triangle area equal to that of a 4.54° equatorial square. His synthesized tides incorporated 
dissipation by means of energy flux across the domain boundary and no tidal loading or self attraction was 
included. 

The numerical ocean tidal models currently available differ by several tens-of-centimeters in many areas 
of the world. The precision altimetric measurements of the sea surface topography now provide the capability 
of independently assessing the accuracy of these models. Simulation studies by Estes (1980) using numerical 
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solutions of the Laplace’s tidal equations to generate the observed measurements (based on simulated 4-day 
SEASAT orbits over a 200-day interval) indicate that when the rms amplitude error is less than 10 cm for 
the M2 tide virtually all the structure is recovered. For rms errors between 12 and 15 cm all large scale 
features are retained, for 15 cm to 20 cm rms errors the high amplitude structure is recovered. These results 
were obtained by utilizing long arc intersection, or crossover point data and simultaneously solving for the 
orbit error parameters with the tide model parameters. 

Mazzega (1985) created a global model of the M2 tide using twenty four days of SEASAT altimetry data, 
the solution was obtained by means of surface spherical harmonics, the results were qualitatively realistic. 
No hydrodynamic equations were used. 

The purpose of this investigation is to develop an objective analysis technique which will serve to ex- 
trapolate tidal amplitudes and phases over entire ocean basins using existing gauge data and the precise 
altimetric measurements which are now beginning to be provided by satellite oceanography. The applicability 
of the extrapolating technique was tested in the Lake Superior basin by Sanchez, Rao and Wolfson (1985). 
Proceeding in stages we now report the results obtained in the Atlantic-Indian oceans using a 6° x 6° grid. 
At this level of resolution it is possible to test the essential features of the method without using real data. 

The method to be used in this investigation requires several distinct steps. First it is necessary to deter- 
mine numerically the stream function and velocity potential orthogonal functions (The Stokes/Helmholtz 
Potentials) which span the space of the basin under consideration. These space functions are then used in 
the Laplace tidal equations to determine the homogeneous solution (normal modes) and the forced solution. 
The latter is obtained by adding the astronomical forcing function modified to include solid-earth tides. 

The velocity potential eigenfunctions obtained as a first step arc used also to extrapolate the surface height 
field over the entire space domain of the given basin and this approach will constitute a distinct and integral 
part of the investigation. 

The theoretical foundation is Proudman’s theory (1918) as formulated by Rao (1966). The theory pro- 
vides the formalism for calculation of the gravitational (first class) normal modes and the rotational (second 
class or Rossby waves) normal models of irregularly shaped basins with realistic bathymetry. 

The method requires the solution of two elliptic partial differential equations with second order operators 
which are simpler than the tidal operator. The boundary conditions correspond to vanishing of the stream 
function and normal derivative of the velocity potential. The elliptic operators are represented numerically 
in finite difference form, the grid used is a Richardson lattice which preserves self-adjointness. The solu- 
tions yield the velocity and surface height fields in terms of orthogonal functions with time-dependent coef- 
ficients. These functions are then substituted into the Laplace’s tidal equations: if the homogeneous equa- 
tions are used one obtains the normal modes; if the forcing terms are included then the forced solution is 
obtained. In both cases the solution is obtained numerically. The surface height field is only dependent on 
the velocity potential orthogonal functions. The expansion coefficients of these functions can be estimated 
in a least-square sense from available selected tidal measurements. 

2. BASIC EQUATIONS 

Free Solutions 


The method of approach was developed originally by Proudman (1918) using a Lagrangian approach. It was reformulated 
by Rao (1966) from the Eulerian point of view. The basic ideas of the method presented below follow Rao’s line of develop- 
ment. The basic equations are the linearized shallow water equations on a rotating plane: 


- f[M] = -g Hh Vq 
4^ + V • M = 0 

at 


( 2 . 1 ) 
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where 


M = HV s (M,N) 


V a (u,v) 
f = 2cosin 0 


h(x,y) = H (x,y)/H. 

H(x,y) is the variable depth of the fluid in equilibrium and H is some constant scaling depth, f is the coriolis parameter, 
V is the horizontal velocity vector, rj the fluctuation of the free surface, g is the apparent gravitational acceleration and V 
is the horizontal gradient operator. []: denotes rotation of the vector through a right angle in the clockwise direction of the 
horizontal plane, i.e. [] = -(k x V), k being a unit vector vertical to the horizontal plane. 

The appropriate boundary conditions to be adjoined to equations (2.1) are 


M ' n = 0 


on the coast, where n is Jhe unit normal to the coast line. 

The transport vector M may be partitioned as follows: 

M = M* + VT 


where 

M* = -h V<)> 

m v = -rvtpi 


( 2 . 2 ) 


(2.3) 


(2.4) 


is the solenoidal part of M while h 'M't' is the irrotational part 


since 


V[h -1 M+] = 0 


VM v = 0 


To complete the determination of <(> and i p, it is necessary to specify the boundary conditions ' n = 0 and ' n 
to ensure that eq. (2.2) is satisfied. In terms of <(> and tp the conditions then are 


= 0 


and 



= 0 


ip = 0 


on the boundaries. 

The divergence of the transport field and the vorticity of the velocity field yield 

V ■ M = — V'hV<)> 


(2.5) 


V ' [h -1 M] = V-h _1 Vrp (2.6) 

♦ 

If M is known as a function of the horizontal coordinates, the left-sides of (2.6) are specified. Then each equation represents 
an inhomogeneous elliptic differential equation with homogeneous boundary conditions given by eq. (2.5), and it is well 
known that such problems possess a unique solution. It is also straight forward to prove that the representation of M as 
given in eqs. (2. 3, 4, 5) is unique. Since M itself is unknown, the procedure then consists of expanding and ip in terms 
of the spectra of the elliptic operators appearing in Eq. (2.6); that is, we seek solutions which satisfy the following equations 
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V • hV<)> y = — A y <(> 

V ' h _1 Vlfly = — HylJJy 


(2.7) 


where <|> y , ip r are the characteristic functions and A y , pi y are the characteristic values associated with the corresponding operators 
V'h V and V‘h -1 V. 

The characteristic functions satisfy the boundary conditions 


,a|r 

an 


h 1 yj y = 0 


= 0 


( 2 . 8 ) 


The condition h _1 ip y = 0 imposes a stronger condition than that required by (2.5). However, the factor h _1 is necessary 
to make the ip y problem self-adjoint. Since the problems (2.7,8) are self-adjoint, the characteristic values A y and p y are real 
and the characteristic functions <(> y and ip y each constitute a complete and internally orthogonal set. Without loss of general- 
ity, the orthogonality condition may be chosen as: 


/h _1 M*'M$dA = Ay/i^pdA = Ac 2 ! 2 ^, 
/h _1 M^ y -M^lA = A y / ip y ippdA = Ac 2 ifd y p, 


(2.9) 


where c 2 = gH and A is the surface area of the basin. 6 y p is the Kronecker delta. We have further defined, in accordance 
with (2.4), 


M+ = — hV<j> y and Mjf = -[Vip y ] 

in the orthogonality condition (2.9). 

The components of the transport vector can now be represented by the sums 


M*=IP M* 

Y y Y’ 

1 * (2.10) 

M v = IQ MS’, 
y r y 

where P y and Q y are the expansion coefficients. The orthogonality conditions ensure that equations (2.10) represent the least 
squares approximations to and M v if and when the sums span the complete spectra of equations (2.7). The height 
field rj is governed by the divergent part of M and the <j> y functions yield a sufficient basis for its representation. A conven- 
ient representation for rj may be taken as: 


9 = 

rj y = c -1 (A y )' /2 <|> y , 


( 2 . 11 ) 


where R y are the expansion coefficients for the rj-field. 

Substitution of equations (2.10) and (2.11) into equations (2.1) and the use of the orthogonality conditions (2.9) yield 
a set of prediction equations for the expansion coefficients: 


^r-lA yP Pp-IB yP Qp-^ y R y = 0 
dt P P 

^ | C yp P p - I D y pQp = 0 

dt p P 


( 2 . 12 ) 
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where the coupling coefficients are given by 


B r p = {M+ y , [M^]} , 
D yp = {M v r , [M^]} , 


A yp = {MP, [M+p]} , 
C y p = {Mjf, [M^]} , 


(2.13) 


The quantity inside the angular brackets represents an inner product of two vectors W and T and is denoted by 


{W,T} = (1/c 2 AH 2 ) / A fh _I WTdA. 

It can be seen from eq. (2.12, 13) and the definition of the inner product that all the coupling coefficients vanish when 
the coriolis parameter f=0 and hence 

d 2 P -> 

— y + = o. 

dt 2 

tf y = (c 2 A y ) ' h , is the non-rotating frequency of oscillation. 

In equations (2.12, 13) the subscript y (or ft) is used as a proxy for a binary index and represents an ordering of the 
characteristic functions <j> y , ip y in some as yet unspecified manner. For convenience we replace the wavenumber index y by 
scalar indices i = 1,2,3 (or j= 1,2,3. . .) and denote 

Ai Ayi, Mj My|, ifj = ifyj 

Pj = Pyj, = Qyi, Rj = Ryj 

A ij s A riffi 

Bij = B yi/Jj (2.14) 

Cij = Cyifo 

D ij = Dy i(Jj 

Using the above notation, we can define column vectors 

P = col (Pj), Q = col (Qj), R = col (Rj) 


and matrices 



A « ( Aij [, B * ( B S | 
C . ( C,j D = { Dj, ) 


<tf> = diagonal 

Equations (2.12) may now be written in the form 

ds + a s = 0 

f\J 

dt 


(2.15a) 


(2.15b) 


(2.16) 


5 



where a is the square matrix 


-A 

-B 

-<tf> 

-C 

-D 

0 

<C> 

0 

0 


(2.17) 


In seeking the solutions for the normal modes, we assume 

s oj e iot 

where o is the normal mode frequency with rotation and i = s — Eq. (2.15) then reduces to 

(oil - ia) s = 0 . (2.18) 

f\J 

In the preceding equation II is the identity matrix, o’s are the characteristic values of the matrix ia. From the definition 
of the coupling coefficients given in (2.13) and the matrix a as defined in (2.17), it is clear that ia exhibits Hermitian 
symmetry and hence the o’s constitute a spectrum of real eigenvalues. In computing the matrix elements in a , the basis 
functions and tp; are chosen to correspond to an ordering of the characteristic values A, and p; arranged so that 
Ai < A 2 < A 3 ... and Mi < M 2 < F 3 • • • • Such an ordering has been chosen since the Aj’s and pi’s have the dimensions of 
(wavenumber) 2 . Hence, at any order of truncation, those 4i’s and ip/s with the largest space scales are taken into account. 

Forced Solution 


In the computation of the forced solution we include the effects of the yielding of the solid earth to tide generating forces. 
The effects due to self-attraction of the tide and tidal loading as well as the frictional effects have not been included. The 
theory allows for the inclusion of these effects but they introduce computational complications and physical uncertainties and 
they are not necessary in the context of the application of the interpolation technique. Equations (2.1) are recast for the forced 
tidal oscillations as 


3M ♦ — 

— - f[M] = — gHhVq', 


where 


dri + 

— + V.M = 0, 

at 


n' = *7 — (1 + k 2 -h 2 )rj. 

-► 

k 2 ,h 2 are the Love numbers and rj is the equilibrium tide height. Let the vector G 
(2.19) can then be written 


(2.19) 


gHh (1 +k 2 — h 2 )Vrj. Equations 


3M ♦ — + 

— - fIM] = -gHhVn + G 

at 

— + V.M = 0, 

at 


( 2 . 20 ) 


Substitution of the expansions (2.10) and (2.11) into equations (2.20) yields, after using (2.13) 

— 1 1 - Z A y flP/j - Z Bwj Qn - R y = (c 2 AH 2 ) /fh _l GM*dA 
dt ft p 


dQ 


- Z CyflP/j - Z D r s Q a = (c 2 AH 2 ) _ 1 /fh _ 1 G'MfflA 
dt p p 

— r + ^y Py = 0 
dt r 7 


( 2 . 21 ) 
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Define the vector F as follows, 


F = 


(^AH 2 ) - 1 / fh - ‘G'MJcIA 
(c 2 AH 2 ) -1 /fh -1 G‘MfdA 


( 2 . 22 ) 


Equations (2.21) can then be written as an inhomogeneous matrix equation 


dS + aS = F, 
dt ~ 

where S and a^ are defined by equations (2.15a) and (2.17). The solution of equation (2.23) is given by 

S = C<e i rt > f„<e "%C _1 F(t) dr, 


(2.23) 


(2.24) 


where 

C is the modal matrix containing the eigenvectors of a arranged in columns and a are the characteristics values of a . <e I0t > 
is a diagonal matrix. 


3. OBJECTIVE ANALYSIS 

The orthogonal functions <f> y form a basis for the expansion of the height field q, as discussed earlier. Since these func- 
tions are characteristic for a particular basin and are described at all points inside the basin (within the resolution of the finite 
difference grid), they can be used as optimal functions for extrapolation of the tidal field over the whole basin, given data 
at some selected points. 

Even through the basis functions are orthogonal, in fitting the data to these functions it is preferable to use linear least- 
squares techniques to determine the coefficients of expansion rather than the orthogonality property since real data are not 
usually available at regularly spaced points. A similar approach was taken in a previous investigation by Rao and Schwab 
(1981) in which they determined the steady circulation in a closed basin for which the appropriate orthogonal functions are 
the ip y functions and Sanchez et al. (1985) used the <|> y functions for tidal extrapolation in Lake Superior. An outline of the 
procedure and the basic equations involved are given below. 

The tidal height field as given by equations (2.11) and (2.15) can be written 

q(x,y,t) = IR /r q y (x,y) cosdt - XR yl q y (x,y) sindt (3.1) 

where R yr and R y j denote the real and imaginary parts of R y . The tidal height field can be expressed also in terms of amplitude 
A(x,y) and phase 0(x,y), that is 


rj(x,y,t) = A(x,y) cos [ot - 0(x,y)] 
Comparing equations (3.1) and (3.2) one obtains the following relations, 

A(x,y) = [(SR yr q y ) 2 + (IR y , r, x ) 2 ] I/2 
0(x,y) = arc tan [-(IR y i q y )/(I R yr q y )] 

Also, from equations (3.1) and (3.2) one obtains the following vector-matrix relations 

[q ] R r = A cos 0 

f\J 

[q ] Rj = A sin 0 


(3.2) 


(3.3) 


(3.4) 
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The components of the vectors A cos 0 and A sin 0 are the “n” available tidal measurements. The matrix [>7 ] is available 
from the solution of the velocity potential eigenfunctions, its “k” columns will correspond to the velocity potential eigenfunc- 
tions chosen to represent the forced solution, its “n” rows correspond to the location of the tidal measurements in the basin 
under consideration. The least squares solution to equations (3.4) is then given by 

Rr = {[n ffo ] r'fo f A 

r\j r\j r\j 

Ri = (fo fin i r'fo f a 

r\j r\j r\j 

T * * 

where [17 ] is the transposed matrix. Having determined a certain number of coefficients R r and Rj, since the functions rj y (x, y) 
are known over the entire basin, one can obtain the amplitude and phase of the tidal height at all points in the basin from 
eq. (3.3). The maximum number of coefficients fe r and Rj that can be determined will be the same as the number of observa- 
tions available. 

4. RESULTS 

Normal Modes and Forced Solutions 

A normal mode solution for the Atlantic and Indian oceans has been obtained by means of a 6° x 6° finite difference 
grid. There are 462 velocity potential and 370 stream function points distributed such as to form a single Richardson lattice. 
A finite difference solution of equations (2.7) subject to the boundary conditions given by equations (2.8) yields the eigen- 
values and eigenvectors for both fields. The eigenfunctions of the 0 solution represent the non-rotating gravitational normal 
modes of the basin, the periods of oscillation of the lowest modes are given in the second column of Table 1. The stream 


TABLE 1. GRAVITATIONAL MODES 

PERIOD (HOURS) ENERGY RATIO 

ROTATIONAL NON-ROTATIONAL (KINETIC/POTENTIAL) 


67.92 

77.68 

1.47 

42.53 

40.98 

0.96 

29.61 

29.00 

0.81 

23.55 

22.17 

0.84 

20.93 

20.73 

0.86 

19.40 

19.20 

1.19 

17.34 

19.08 

1.03 

15.73 

16.52 

1.11 

14.44 

14.56 

0.86 

13.93 

14.17 

1.03 

12.99 

13.42 

1.07 

12.55 

12.94 

1.27 

11.99 

12.74 

1.09 

11.13 

11.72 

1.14 

10.93 

11.47 

1.14 

10.62 

11.12 

1.15 

10.39 

10.60 

1.21 

9.97 

10.26 

1.16 

9.52 

9.99 

1.19 

9.21 

9.26 

1.00 


cos 0, 
sin 0, 


(3.5) 
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function modes generate the vorticity component of the flow field in the rotating case for the gravitational and rotational species, 
they are more dominant for the quasi-geostrophic rotational modes than for the gravitational modes. The normal mode solu- 
tion was obtained by including the lowest 150 eigenfunctions from each field (<f> and ip ) into the dynamical equations, the 
solution of the eigenvalue problem posed by equation (2.18) yields the normal modes. The normal modes fall into two distinct 
categories: the inertia-gravitational modified by rotation and the rotational modified by divergence. If the secular determinant 
is truncated at a size n x n, then one obtains 2n/3 gravitational modes and n/3 rotational modes, in pairs (plus and minus 
the same value). In our case there are 300 gravitational modes and 150 rotational modes. The rotational modes are characterized 
by high ratios of kinetic to potential energy, these modes will vanish in a basin of constant depth and Coriolis parameter. 
The rotational normal modes are extremely sensitive to the resolution of the shape of the basin and the bottom topography 
and the convergence of their roots is more complicated than for gravitational modes. A detailed discussion of the nature of 
the two types of modes is given by Platzman (1975) and Rao and Schwab (1976). The gravitational modes are by far the 
most important in the context of the diumal and semidiurnal tides as shown by Platzman (1984). 

Table 1 gives a list of the first 20 gravitational modes in order of decreasing period, for each mode the following quantities 
are given: the rotational and non-rotational periods in hours, the ratio of energies (kinetic/potential). 


Table 2 lists the 20 most powerful modes for the M2 and K1 tidal components in order of decreasing power. The entries 
in the table give the power ranking, the rotational period in hours, the percentage of the total power for that mode and the 
percentage of the total power contributed by all the modes included up to that point. 


TABLE 2. POWER SPECTRUM FOR THE M2 AND K1 COMPONENTS 

M2 K1 

RANK PERIOD POWER (%) PERIOD POWER (%) 



(HOURS) 

EACH 

SUM 

(HOURS) 

EACH 

SUM 

1 

12.99 

7.67 

7.67 

23.55 

12.79 

12.79 

2 

10.62 

7.33 

15.01 

13.93 

8.70 

21.49 

3 

13.93 

7.31 

22.32 

17.34 

8.19 

29.68 

4 

11.99 

6.00 

28.33 

20.93 

8.01 

37.70 

5 

14.44 

5.34 

33.67 

29.61 

3.52 

41.23 

6 

10.39 

3.84 

37.51 

15.73 

3.47 

44.70 

7 

12.55 

3.61 

41.13 

19.40 

3.37 

48.08 

8 

9.21 

3.44 

44.58 

12.99 

2.92 

51.00 

9 

19.40 

2.95 

47.53 

10.93 

2.82 

53.82 

10 

17.34 

2.74 

50.28 

11.13 

2.62 

56.45 

11 

9.52 

2.68 

52.96 

9.21 

2.15 

58.60 

12 

8.92 

2.62 

55.59 

67.92 

2.04 

60.65 

13 

11.13 

2.52 

58.12 

14.44 

1.62 

62.28 

14 

8.30 

2.41 

60.53 

8.80 

1.61 

63.90 

15 

9.97 

1.62 

62.16 

8.30 

1.58 

65.49 

16 

8.01 

1.55 

63.71 

11.99 

1.46 

66.95 

17 

23.55 

1.54 

65.25 

9.97 

1.36 

68.31 

18 

8.70 

1.44 

66.70 

12.55 

1.34 

69.66 

19 

8.09 

1.43 

68.14 

10.39 

1.34 

71.00 

20 

7.82 

1.40 

69.54 

8.40 

1.14 

72.15 


The three most powerful modes in the M2 spectrum are the modes with periods of 12.99 hours (7.67%), 10.62 hours 
(7.33%)and 13.93 hours (7.31%). FortheKl spectrum the three most powerful modes have periods of 23.55 hours (12.79%), 
13.93 hours (8.70%) and 17.34 hours (8.19%). Figures 1 to 5 below give the structure (amplitude and phase) for these 5 
modes. The amplitudes have been normalized to a maximum value of 100, the amplitude contours are shown by the solid 
lines and the dotted lines correspond to the phases. 
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Figure 3. Normal Mode, Period = 



Figure 4. Normal Mode, Period 



12.99 Hours 







Figure 5. Normal Mode, Period = 10.62 Hours 


Figures 6 and 7 show the M2 and K1 tidal solutions. These solutions were obtained by including the first 150 eigenfunc- 
tions from the and ip solutions. The forced solution for rj then contains 150 coefficients R y . The contours of equal amplitude 
and phase are given by the solid and dotted lines respectively. The arrows indicate the sense of progression for high and 
low water. 

The solution for the M2 tide shows an amphidromic point south of Greenland and a double amphidromic system in the 
central North Atlantic. There is an amphidrome in the South Atlantic and several amphidromes in the periphery of the Indian 
Ocean, the central Indian Ocean shows the usual anti-amphidromic high. The K1 tide shows a simpler structure with am- 
phidromic points in the central North and South Atlantic and two amphidromes in the Indian Ocean. The M2 solution has 
an amplitude rms of 51.6 cm, the corresponding value for K1 is 26 cm. 

Figures 8 and 9 show the approximate location of the amphidromes for the M2 and K1 tidal solutions as obtained by 
different investigators. The different models agree fairly well with each other in certain areas and show wide variations in 
other regions. The results of our model are reasonable, undoubtedly much of the variations can be attributed to poor resolution 
due to the coarseness of the 6° x 6° grid as well as the absence of frictional effects, self-attraction and tidal loading. 
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Figure 8. Location of Amphidromes for the M2 Tide 

o: Schwiderski (1979) a: Accad and Pckeris (1978) 

O: Parke and Hendershott (1980) □ : This Investigation 



Figure 9. Location of Amphidromes for the K1 Tide 

o: Schwiderski (1981) o: Parke and Hendershott (1980) □ : This Investigation 
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Objective Analysis Evaluation 


In the objective analysis of the tide we attempt to use the data wherever they are available and fit co-range and co-tidal 
lines to describe the tide over the entire basin using the velocity potential eigenfunctions in a least squares sense. Since the 
actual tidal distribution over the entire basin is unknown, it is not possible to answer the question of how well this extrapola- 
tion procedure has worked. We have chosen to use the theoretically simulated tidal field as the “true solution” and compare 
the resulting tidal pattern obtained from a selected set of data points. The theoretical simulation also serves as a guide in 
choosing the most energetic modes to be used in the objective analysis. Ideally, in the least squares procedure it is preferable 
to use fewer functions than there are available data points to ensure a smooth fit and employ an a priori knowledge of the 
most dominant modes in making this selection. 

In order to carry out the testing of the objective analysis two sets of data points have been created. In one case every 
fourth point in the <)> field was defined as a data point, yielding the 116 data points shown in Figure 10. In the second case 
every sixth point was taken, to obtain the 77 data points shown in Figure 1 1 . The tidal amplitudes and phases at these points 
were extracted from the theoretical simulation of the M2 tide using 150 modes. 

A comparison of the objectively analyzed amplitudes and phases with the “true” values at all grid points is shown in 
Tables 3 and 4. The tables list the percentage of total grid points (462) where the analyzed and true values of amplitude 
and phase agree with each other within the limits indicated for the different cases corresponding to the number of modes 
chosen in the analysis. Also given in the table are the rms errors of the amplitude over the entire basin and only at the data points. 

In addition the two data sets were objectively analyzed using spherical harmonics. The results are shown in Tables 5 and 
6. These tables list the same quantities as in Tables 3 and 4 for the different cases corresponding to the degree of the expansion 
chosen for the fit. 
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Figure 11. Data Set (77 Points) 


TABLE 3. COMPARISON BETWEEN OBJECTIVELY ANALYZED AND THEORETICAL FIELDS 

.M2 COMPONENT. (116 data points) 

Coefficients Amplitude (%) Phase (%) RMS (cm) 

(% of Total Power) <5% <10% <20% <5° <10° <20° TOTAL DATA 


14 

(60.5) 

14 

27 

48 

29 

48 

71 

11.6 

9.4 

29 

(78.5) 

35 

55 

79 

51 

75 

90 

5.4 

3.6 

39 

(84.4) 

37 

60 

83 

56 

78 

94 

7.1 

2.7 

48 

(88.0) 

47 

69 

88 

65 

84 

95 

7.0 

2.0 

58 

(90.8) 

53 

78 

92 

68 

87 

97 

4.2 

1.6 

68 

(92.8) 

56 

74 

89 

73 

84 

92 

9.1 

0.9 

80 

(94.7) 

58 

73 

88 

76 

89 

96 

19.5 

0.6 
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TABLE 4. COMPARISON BETWEEN OBJECTIVELY ANALYZED AND THEORETICAL FIELDS 

.M2 COMPONENT. (77 data points) 

Coefficients 
(% of Total Power) 

Amplitude (%) 

<5% <10% <20% 

A 

U\ 

o 

Phase (%) 
<10° 

<20° 

RMS (cm) 
TOTAL DATA 

8 

(44.5) 

8 

18 

31 

9 

18 

37 

18.8 

17.2 

14 

(60.5) 

14 

27 

51 

30 

50 

70 

13.6 

9.8 

19 

(68.1) 

16 

34 

58 

35 

55 

79 

12.4 

7.8 

24 

(74.0) 

26 

43 

68 

36 

60 

83 

12.9 

4.7 

29 

(78.5) 

35 

55 

78 

44 

63 

86 

10.4 

2.9 

32 

(80.4) 

32 

57 

78 

43 

63 

86 

8.6 

2.5 

34 

(81.6) 

34 

55 

78 

46 

68 

88 

15.5 

1.9 

39 

(84.4) 

32 

52 

76 

46 

67 

86 

35.5 

1.7 


TABLE 5. OBJECTIVE ANALYSIS USING SPHERICAL HARMONICS .M2 COMPONENT. 

(116 data points) 


Degree 

(No. of Coeff.) 


Amplitude (%) 

<5% <10% <20% 


Phase (%) 

<5° <10° <20° 


RMS (cm) 
TOTAL DATA 


(18) 

5 

11 

23 

5 

11 

28 

30.0 

33.2 

(32) 

6 

13 

28 

8 

18 

36 

26.4 

29.4 

(50) 

7 

15 

28 

15 

30 

52 

21.8 

19.5 

(72) 

13 

24 

47 

24 

44 

66 

19.2 

11.6 

(98) 

15 

30 

56 

27 

51 

71 

21.7 

7.8 


TABLE 6. OBJECTIVE ANALYSIS USING SPHERICAL HARMONICS. M2 COMPONENT. 

(77 data points) 


Degree 

(no. of coeff.) 


Amplitude (%) 

<5% <10% <20% 


Phase (%) 

<5° <10° <20° 


RMS (cm) 
TOTAL DATA 


(18) 

9 

14 

27 

3 

9 

18 

25.2 

21.5 

(32) 

10 

18 

39 

14 

26 

42 

21.9 

15.0 

(50) 

11 

22 

37 

17 

30 

54 

19.2 

13.2 

(72) 

16 

28 

50 

21 

35 

51 

34.5 

7.9 








Tables 7 and 8 show the results obtained for the objective analysis of the K1 tide. When using the data set with 116 points 
the best fit to the M2 tide is obtained when estimating the coefficients corresponding to the 58 most powerful modes, this 
case yields a 4.2 cm RMS over the entire 462 point 4 field. For the K1 tide the best fit corresponds to the case when 35 
coefficients are estimated wnich yields a total RMS of 1.9 cm. For the data set containing 77 points the best fit to the M2 
tide is obtained when estimating the coefficients corresponding to the 32 most powerful modes which case gives a total RMS 
of 8.6 cm. For this data set the best fit to the K1 tide yields a 3.3 cm total RMS by estimating 34 coefficients. The number 
of coefficients which yields the best fit is a function of the number of data points and their distribution in space in relation 
to the structure of the particular tidal component. Therefore, given a certain data set in the domain of an ocean basin, the 
theoretical simulation will help to choose the optimum number of functions to use in the objective analysis. Also by using 
the results of the objective analysis of the theoretical solution it is possible to determine the areas where additional data 


TABLE 7. COMPARISON BETWEEN OBJECTIVELY ANALYZED AND THEORETICAL FIELDS. 

K1 COMPONENT. (116 data points) 

Coefficients 
(% of Total Power) 

Amplitude (%) 

<5% <10% <20% 

A 

U\ 

o 

Phase (%) 
<10° 

<20° 

RMS (cm) 
TOTAL DATA 

20 

(72.1) 

29 

48 

68 

43 

71 

94 

3.4 

3.2 

35 

(83.0) 

44 

69 

89 

69 

89 

97 

1.9 

1.2 

40 

(85.1) 

51 

73 

89 

73 

90 

97 

2.1 

1.0 

50 

(88.3) 

62 

81 

93 

73 

93 

99 

2.2 

0.8 

70 

(92.7) 

63 

84 

94 

78 

92 

97 

2.2 

0.3 

100 

(96.9) 

69 

81 

94 

79 

88 

96 

8.1 

0.07 


TABLE 8. COMPARISON BETWEEN OBJECTIVELY ANALYZED AND THEORETICAL FIELDS. 

K1 COMPONENT. (77 data points) 

Coefficients 
(% of Total Power) 

Amplitude (%) 

<5% <10% <20% 

o 

V 

Phase (%) 
<10° 

<20° 

RMS (cm) 
TOTAL DATA 

20 

(72.1) 

24 

42 

66 

42 

66 

90 

4.4 

3.7 

23 

(75.3) 

26 

45 

70 

44 

73 

89 

3.6 

3.5 

32 

(81.6) 

40 

61 

81 

54 

78 

94 

5.9 

1.1 

34 

(82.6) 

39 

62 

81 

53 

79 

94 

3.3 

1.0 

36 

(83.5) 

40 

66 

83 

57 

79 

96 

3.5 

0.8 

45 

(86.8) 

50 

72 

89 

66 

83 

95 

6.8 

0.5 
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measurements are most needed as shown in Figures 12 and 13. Figure 12 shows the absolute values of the differences in 
amplitudes between the theoretical M2 tide using 150 coefficients and the best fit using 58 coefficients and 116 data points. 
Figure 13 shows the corresponding values for the K1 tide when using a best fit of 35 coefficients. The numbers are rounded 
to the nearest digit to facilitate display. In the M2 case there are areas in the north-east Atlantic with higher than average 
RMS. The highest discrepancies occur at some isolated points off the coast of Labrador, Greenland, northern Australia and 
east Antarctica. The K1 fit shows some high values in the north-east Atlantic and the Patagonian shelf as well as in the Antarc- 
tic area in the Indian Ocean. 

The results of the objective analysis of the M2 tide using spherical harmonics given in Tables 5 and 6 indicate a best 
fit when using an expansion of degree 5 with the 116 point data set and a best fit of degree 4 for the 77 point data set. Both 
cases yielded a total RMS of 19.2 cm, well above the values of 4.2 cm and 8.6 cm obtained when using the <|> eigenfunctions. 
The K1 tide was analyzed also using spherical harmonics, the results given in Tables 9 and 10 show a best fit of 3.3 cm 
when using an expansion of degree 5 with the 1 16 point data set and a best fit of 5.9 cm when using an expansion of degree 
4 with the 77 point data set. 

Table 1 1 contains a summary and comparison between the best cases using both spherical harmonics and the eigenfunc- 
tions of the velocity potential which are referred to as Proudman functions. The Proudman functions yield better results uniformly 
when comparing corresponding tidal components and data sets, the improvement is more pronounced when fitting the M2 
constituent which has a more complex structure. We shall point out that we have not tested every possible combination of 
modes when estimating the coefficients to fit the tide in the different cases. An exhaustive search might yield still better results 
for the Proudman functions. 



Figure 12. M2 Tide: Amplitude Differences (cm) Between Theoretical Solution Using 150 Coefficients and 

58 Coefficient Solution from 116 Point Data Set. 
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Figure 13. K1 Tide: Amplitude Differences (cm) Between Theoretical Solution Using 150 Coefficients and 

35 Coefficient Solution from 116 Point Data Set. 


TABLE 9. OBJECTIVE ANALYSIS USING SPHERICAL HARMONICS. 
K1 COMPONENT. (116 data points). 


Degree 

(No. of Coeff.) 



TABLE 10. OBJECTIVE ANALYSIS USING SPHERICAL HARMONICS. 
K1 COMPONENT. (77 data points). 


Degree 

(No. of Coeff.) 



Amplitude (%) 

< 5 % < 10 % < 20 % 


Phase (%) 

< 5 ° < 10 ° < 20 ° 


RMS (cm) 
TOTAL DATA 


7 

16 

33 

9 

19 

37 

10.8 

11.4 

13 

29 

52 

22 

43 

70 

6.2 

5.6 

19 

43 

70 

39 

66 

90 

5.9 

2.6 

24 

45 

70 

39 

63 

87 

10.2 

1.8 
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TABLE 11. COMPARISON BETWEEN PROUDMAN FUNCTIONS AND SPHERICAL HARMONICS 

Amplitude RMS 
(cm) (% of mean value of 

Tidal Component Data Set Method theoretical solution) 


M2 

116 

Proudman functions 

4.2 

8% 

M2 

116 

Spherical harmonics 

19.2 

37% 

M2 

77 

Proudman functions 

8.6 

17% 

M2 

77 

Spherical harmonics 

19.2 

37% 

K1 

116 

Proudman functions 

1.9 

7% 

K1 

116 

Spherical harmonics 

3.3 

13% 

K1 

77 

Proudman functions 

3.3 

13% 

K1 

77 

Spherical harmonics 

5.9 

23% 


5. CONCLUSIONS 

A method has been developed to objectively analyze tidal data to map the tidal components over an entire basin. The 
technique has been demonstrated in an application to the Lake Superior basin by Sanchez, Rao and Wolfson (1985). Now 
it has been applied to the Atlantic-Indian ocean basin using a 6° x 6° finite difference grid. 

The computed normal modes agree fairly well with those obtained by Platzman (1975) using a different approach. The 
forced solutions for the M2 and K1 tidal constituents were computed also, these theoretical simulations show the same general 
pattern obtained by other modelers, i.e. Platzman (1984), Schwiderski (1983), Parke and Hendershott (1980) and Accad and 
Pekeris (1978). The theoretical solutions were obtained by including 150 eigenfunctions of the velocity potential, the resulting 
height and phase fields were used as data to test the objective analysis technique and as a guide in choosing the most energetic 
modes. 

The results of the objective analysis show that it is possible to recover the M2 and K1 components with a degree of ac- 
curacy well within the error bounds of present day satellite techniques and with amplitude errors smaller than the differences 
existing between different numerical models in many areas of the world oceans. The Proudman functions performed con- 
sistently better than spherical harmonics in the recovery of the tides especially for the M2 constituent with its complex structure. 

Future applications of the method should consider its extension to the Pacific basin, the adoption of a grid with finer resolution, 
the incorporation of real data from gauges and satellite measurements and the simulation and analysis of other tidal consti- 
tuents. The method is especially suited for the latter since the eigenfunctions of the velocity potential and the normal modes 
of a particular basin are dependent only on its shape and depth and they have to be computed only once. 
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